#Taylor N. Carlson
#1.27.16
#Political Chameleons Study 2 Replication Code

#This file includes the replication code for analyses in Study 2 of "Political Chameleons: An Exploration of Conformity in Political Discussions"
#The file first shows the analyses for summary statistics, manipulation checks, and treatment main effects. Second, it includes the code that generated the figures in the paper. Third, it includes the code for the regression models

rm(list=ls())

#Load packages 
library(stargazer)
library(xtable)

#Set working directory
#setwd("/Users/Taylor/Box Sync/Contentious Social Interaction Projects/Political Chameleons/Final Manuscript Materials/Polcham Replication Materials/")

#Load in data
pol<-read.csv("study2data.csv")
dim(pol)
pol2<-subset(pol, !is.na(pol$pressure)) #Just those who took the posttest
dim(pol)
dim(pol2)

#Summary Statistics
summary(pol$conformcount)
summary(pol2$conformcount_strict)


#Manipulation Checks

     #Faux vs. critical questions

t.test(pol$fauxchange, pol$criticalchange)
t.test(pol2$fauxchange, pol2$criticalchange)
     #Significantly different
     
#Checking Differences Between Randomly Assigned Question Orders

#Conformity
	#Potential
		t.test(pol$conformcount ~ pol$order)
		t.test(pol2$conformcount ~ pol2$order)

	#Pure
		t.test(pol2$conformcount_strict ~ pol2$order)


###Differences Between Conditions

#Conformity 
     #Potential
     t.test(pol$conformcount ~ pol$condition)
     
     #Do a permutation test for the potential conformity condition, given that permutation tests are sometimes more appropriate in small samples
     #We have 63 observations, 34 in the treatment group, 29 in the control group. 
     #Create a dataframe where one column is the label (treatment) and one column is the data (number of times the participant conformed, by pure standards)
     dat1<-pol[, c("condition", "conformcount")]
     head(dat1)
     dim(dat1)
     
     set.seed(1414)
     onePermute<-function(dat1){
     	#reshuffle the labels
     	labels<-c(rep(0, 29), rep(1, 34))
     	labels<-labels[sample(1:63,63)]
     	#replace them in the data
     	dat1[,"condition"]<-labels
     	#aggregate
     	means<-aggregate(conformcount~condition, FUN=mean, data=dat1)
     	meanDiff<-means[2,2]-means[1,2] #treatment mean-control mean
     }

     #Now do this 10,000 times to create a null distribution
     meanDifference_potential<-NA
     for(i in 1:10000){
     	meanDifference_potential[i]<-onePermute(dat1)
     }
     
     #Plot histogram
     hist(meanDifference_potential, main="Histogram of Permutation Test Results\nPotential Conformity (10,000 Simulations)", xlab="Difference of Means (Treatment-Control)", ylab="Frequency", breaks=50, ylim=c(0,1500), xlim=c(-1,1))
     #Plot a line showing what we got in the actual data
     abline(v=mean(pol$conformcount[pol$condition==1])-mean(pol$conformcount[pol$condition==0]), col="red")
     
     #Calculate p-value
mean(meanDifference_potential>mean(pol$conformcount[pol$condition==1])-mean(pol$conformcount[pol$condition==0]))
		#p=.001.  This is the proportion of permuted differences greater than the observed difference.

     #Pure
     t.test(pol2$conformcount_strict ~ pol2$condition)
     	#p=.1047
     	
     #Do a permutation test for the pure conformity condition, given the smaller sample size
     
     #We have 46 observations, 23 in the treatment group, 23 in the control group. 
     #Create a dataframe where one column is the label (treatment) and one column is the data (number of times the participant conformed, by pure standards)
     dat<-pol2[, c("condition", "conformcount_strict")]
     head(dat)
     dim(dat)
     
     set.seed(1515)
     onePermute<-function(dat){
     	#reshuffle the labels
     	labels<-c(rep(0, 23), rep(1, 23))
     	labels<-labels[sample(1:46,46)]
     	#replace them in the data
     	dat[,"condition"]<-labels
     	#aggregate
     	means<-aggregate(conformcount_strict~condition, FUN=mean, data=dat)
     	meanDiff<-means[2,2]-means[1,2] #treatment mean-control mean
     }

     #Now do this 10,000 times to create a null distribution
     meanDifference<-NA
     for(i in 1:10000){
     	meanDifference[i]<-onePermute(dat)
     }
     
     #Plot histogram
     hist(meanDifference, main="Histogram of Permutation Test Results\nPure Conformity (10,000 Simulations)", xlab="Difference of Means (Treatment-Control)", ylab="Frequency", breaks=50, ylim=c(0,1500), xlim=c(-1,1))
     #Plot a line showing what we got in the actual data
     abline(v=mean(pol2$conformcount_strict[pol2$condition==1])-mean(pol2$conformcount_strict[pol2$condition==0]), col="red")
     
     #Calculate p-value
mean(meanDifference>mean(pol2$conformcount_strict[pol2$condition==1])-mean(pol2$conformcount_strict[pol2$condition==0]))
		#p=.0302.  This is the proportion of permuted differences greater than the observed difference. 
     

###Check for differences based on those who might’ve known what was going on in the study

	#for pol$remove, 1=maybe knew too much; 0=no idea what was going on
	t.test(pol$conformcount ~ pol$remove)
	t.test(pol2$conformcount_strict ~ pol2$remove)


#See if results hold excluding those who knew too much
t.test(pol$conformcount[pol$remove==0]~ pol$condition[pol$remove==0])
t.test(pol$conformcount_strict[pol$remove==0]~ pol$condition[pol$remove==0])
t.test(pol2$conformcount_strict[pol2$remove==0]~ pol2$condition[pol2$remove==0])

#Check for results based on the actors
chisq.test(pol$actor1, pol$conformcount, simulate.p.value=TRUE)
chisq.test(pol$actor2, pol$conformcount, simulate.p.value=TRUE)
chisq.test(pol2$actor1, pol2$conformcount_strict, simulate.p.value=TRUE)
chisq.test(pol2$actor2, pol2$conformcount_strict, simulate.p.value=TRUE)
	#Note that we use simulated p values because the sample size is too small for chi square tests and it otherwise generates warning messages

#Are any of these results conditional on demographics?

#Gender:
gender_conform<-aov(pol$conformcount ~ pol$female*pol$condition)
summary(gender_conform) #only a main effect of condition, no interaction or main effect of gender

gender_conform_pure<-aov(pol2$conformcount_strict ~ pol2$female*pol2$condition)
summary(gender_conform_pure) #only a main effect of condition, no interaction or main effect of gender

#Race
race_conform<-aov(pol$conformcount ~ pol$white*pol$condition)
summary(race_conform) #only a main effect of condition

race_conform_pure<-aov(pol2$conformcount_strict ~ pol2$white*pol2$condition)
summary(race_conform_pure) 

summary(lm(conformcount ~ condition + condition*polknow, data=pol)) #no effect of knowledge
summary(lm(conformcount ~ condition + condition*engage2012, data=pol)) #no effect of engage2012
summary(lm(conformcount ~ condition + condition*payattn, data=pol)) #no effect of pay attention
summary(lm(conformcount ~ condition + condition*attach, data=pol)) #no effect of anything
summary(lm(conformcount ~ condition + condition*ideology, data=pol)) #no effect of anything 

#-----Figures-----#

#Standard Conformity by Condition 

conformplot<-c(mean(pol$conformcount[pol$condition==0]), mean(pol$conformcount[pol$condition==1]))

#Plot this vector

barplot(conformplot, names.arg=c("Control", "Treatment"), col=c("lightskyblue", "white"), main="Potential Conformity by Condition", ylab="Mean Frequency of Conformity", xlab="Condition", ylim=c(0,3))

#Calculate 2x standard error of the mean

conSE95<-2*(sd(pol$conformcount[which(pol$condition==0)])/sqrt(length(pol$conformcount[which(pol$condition==0)])))

treatSE95<-2*(sd(pol$conformcount[which(pol$condition==1)])/sqrt(length(pol$conformcount[which(pol$condition==1)])))

#make a vector with mean, mean+SE, mean-SE for control and treatment groups

con95<-c(mean(pol$conformcount[which(pol$condition==0)]), mean(pol$conformcount[which(pol$condition==0)])+conSE95, mean(pol$conformcount[which(pol$condition==0)])-conSE95)

treat95<-c(mean(pol$conformcount[which(pol$condition==1)]), mean(pol$conformcount[which(pol$condition==1)])+treatSE95, mean(pol$conformcount[which(pol$condition==1)])-treatSE95)

#use cbind(con95, treat95) for the table in the rest of the code

cbind(con95, treat95)

#Add error bars
#make a table with the upper bound of the CI as the second row, lower bound as third row

conditiontable95<-cbind(con95, treat95)

h<-conditiontable95[2,] #defining the high bound
l<-conditiontable95[3,] #defining the low bound
x<-seq(1:ncol(conditiontable95))
ind<-c(0.7, 1.9, 3.1, 4.3) #finding center points of each bar

for(i in 1:ncol(conditiontable95)){
	lines(c(ind[i],ind[i]), c(l[i], h[i]), lwd=1, lty=1, col=1)
	}
	
#Strict Conformity by Condition 

conformplot2<-c(mean(pol2$conformcount_strict[pol2$condition==0]), mean(pol2$conformcount_strict[pol2$condition==1]))

#Plot this vector

barplot(conformplot2, names.arg=c("Control", "Treatment"), col=c("lightskyblue", "white"), main="Pure Conformity by Condition", ylab="Mean Frequency of Conformity", xlab="Condition", ylim=c(0,3))

#Calculate 2x standard error of the mean

conSE95<-2*(sd(pol2$conformcount_strict[which(pol2$condition==0)])/sqrt(length(pol2$conformcount_strict[which(pol2$condition==0)])))

treatSE95<-2*(sd(pol2$conformcount_strict[which(pol2$condition==1)])/sqrt(length(pol2$conformcount_strict[which(pol2$condition==1)])))

#make a vector with mean, mean+SE, mean-SE for control and treatment groups

con95<-c(mean(pol2$conformcount_strict[which(pol2$condition==0)]), mean(pol2$conformcount_strict[which(pol2$condition==0)])+conSE95, mean(pol2$conformcount_strict[which(pol2$condition==0)])-conSE95)

treat95<-c(mean(pol2$conformcount_strict[which(pol2$condition==1)]), mean(pol2$conformcount_strict[which(pol2$condition==1)])+treatSE95, mean(pol2$conformcount_strict[which(pol2$condition==1)])-treatSE95)

#use cbind(con95, treat95) for the table in the rest of the code

cbind(con95, treat95)

#Add error bars
#make a table with the upper bound of the CI as the second row, lower bound as third row

conditiontable95<-cbind(con95, treat95)

h<-conditiontable95[2,] #defining the high bound
l<-conditiontable95[3,] #defining the low bound
x<-seq(1:ncol(conditiontable95))
ind<-c(0.7, 1.9, 3.1, 4.3) #finding center points of each bar

for(i in 1:ncol(conditiontable95)){
	lines(c(ind[i],ind[i]), c(l[i], h[i]), lwd=1, lty=1, col=1)
	}

#Create a cumulative growth plot to show the number of people conforming at least once (potential)

#First, total conform count 
x<-table(pol$conformcount)
#remove 0s (first column)
x<-x[-1]
#Calculate the cumulative sum for and turn into a proportion
conformcount_total<-cumsum(x)/length(pol$conformcount)
conformcount_total2<-c(0,conformcount_total, rep(conformcount_total[4], 6))
conformcount_total2

#Now, just conform count for treatment condition
x.treat<-table(pol$conformcount[pol$condition==1]) 
x.treat<-x.treat[-1]
x.treat
#Calculate the cumulative sum for and turn into a proportion
conformcount_treat<-cumsum(x.treat)/length(pol$conformcount[pol$condition==1])
conformcount_treat2<-c(0,conformcount_treat, rep(conformcount_treat[4], 6))
conformcount_treat2

#Now, just conform count for control condition
x.control<-table(pol$conformcount[pol$condition==0]) 
x.control<-x.control[-1]
x.control
#Calculate the cumulative sum for and turn into a proportion
conformcount_control<-c(cumsum(x.control)/length(pol$conformcount[pol$condition==0]))
conformcount_control2<-c(0,conformcount_control, rep(conformcount_control[3], 7))
conformcount_control2

#Total conforming at each amount
conformcount_raw<-table(pol$conformcount)
conformcount_raw<-c(conformcount_raw, rep(0,6))
conformcount_raw<-conformcount_raw/length(pol$conformcount)
conformcount_raw

#Total conforming at each amount (treatment)
conformcount_raw_treat<-table(pol$conformcount[pol$condition==1])
conformcount_raw_treat<-c(conformcount_raw_treat, rep(0,6))
conformcount_raw_treat<-c(conformcount_raw_treat/length(pol$conformcount[pol$condition==1]))
conformcount_raw_treat

#Total conforming at each amount (control)
conformcount_raw_control<-table(pol$conformcount[pol$condition==0])
conformcount_raw_control<-c(conformcount_raw_control, rep(0,7))
conformcount_raw_control<-c(conformcount_raw_control/length(pol$conformcount[pol$condition==0]))
conformcount_raw_control



#Now plot
xaxis<-c(0:10)
par(mar=c(5,4,4,4)+.1)
plot(xaxis, conformcount_total2, type="l", col="slategray4", lwd=1, ylim=c(0,1),xaxt="n", main="Frequency of Potential Conformity", xlab="Total Number of Times Participant Conformed", ylab="Proportion of Participants")
lines(xaxis, conformcount_treat2, col="steelblue1", lwd=1, lty=5)
lines(xaxis, conformcount_control2, col="black", lwd=1, lty=3)
lines(xaxis, conformcount_raw, col="slategray4", lwd=3.5, lty=1)
lines(xaxis, conformcount_raw_treat, col="steelblue1", lwd=3.5, lty=5)
lines(xaxis, conformcount_raw_control, col="black", lwd=3.5, lty=3)
axis(1, at=c(0,1,2,3,4,5,6,7,8,9,10), labels=xaxis)
legend(0,1, legend=c("All Participants", "Treatment", "Control"), lty=c(1, 5,3), lwd=rep(1,3), col=c("slategray4", "steelblue1", "black"), title="Cumulative Distribution", bty="n")
legend(6.5,.3, legend=c("All Participants", "Treatment", "Control"), lty=c(1,5,3), lwd=c(3.5, 2.5, 3.5), col=c("slategray4", "steelblue1", "black"), title="Raw Distribution", bty="n")

#Create a cumulative growth plot to show the number of people conforming at least once (pure)

#First, total conform count 
x<-table(pol2$conformcount_strict)
#remove 0s (first column)
x<-x[-1]
x
#Calculate the cumulative sum for and turn into a proportion
conformcount_total_strict<-cumsum(x)/length(pol2$conformcount_strict)
conformcount_total_strict2<-c(0,conformcount_total_strict, rep(conformcount_total_strict[2], 8))
conformcount_total_strict2


#Now, just conform count for treatment condition
x.treat.s<-table(pol2$conformcount_strict[pol2$condition==1]) 
x.treat.s<-x.treat.s[-1]
x.treat.s
#Calculate the cumulative sum for and turn into a proportion
conformcount_treat_strict<-cumsum(x.treat.s)/length(pol2$conformcount_strict[pol2$condition==1])
conformcount_treat_strict2<-c(0,conformcount_treat_strict, rep(conformcount_treat_strict[2], 8))
conformcount_treat_strict2

#Now, just conform count for control condition
x.control.s<-table(pol2$conformcount_strict[pol2$condition==0]) 
x.control.s<-x.control.s[-1]
x.control.s
#Calculate the cumulative sum for and turn into a proportion
conformcount_control_strict<-c(cumsum(x.control.s)/length(pol2$conformcount_strict[pol2$condition==0]))
conformcount_control_strict2<-c(0,conformcount_control_strict, rep(conformcount_control_strict[2], 8))
conformcount_control_strict2


#Total conforming at each amount
conformcount_raw_strict<-table(pol2$conformcount_strict)
conformcount_raw_strict<-c(conformcount_raw_strict, rep(0,8))
conformcount_raw_strict<-conformcount_raw_strict/length(pol$conformcount_strict)
conformcount_raw_strict

#Total conforming at each amount (treatment)
conformcount_raw_treat_strict<-table(pol2$conformcount_strict[pol2$condition==1])
conformcount_raw_treat_strict<-c(conformcount_raw_treat_strict, rep(0,8))
conformcount_raw_treat_strict<-c(conformcount_raw_treat_strict/length(pol2$conformcount_strict[pol2$condition==1]))
conformcount_raw_treat_strict

#Total conforming at each amount (control)
conformcount_raw_control_strict<-table(pol2$conformcount_strict[pol2$condition==0])
conformcount_raw_control_strict<-c(conformcount_raw_control_strict, rep(0,8))
conformcount_raw_control_strict<-c(conformcount_raw_control_strict/length(pol2$conformcount_strict[pol2$condition==1]))
conformcount_raw_control_strict


#Now plot
xaxis<-c(0:10)
par(mar=c(5,4,4,4)+.1)
plot(xaxis, conformcount_total_strict2, type="l", col="slategray4", lwd=1, ylim=c(0,1),xaxt="n", main="Frequency of Pure Conformity", xlab="Total Number of Times Participant Conformed", ylab="Proportion of Participants")
lines(xaxis, conformcount_treat_strict2, col="steelblue1", lwd=1, lty=5)
lines(xaxis, conformcount_control_strict2, col="black", lwd=1, lty=3)
lines(xaxis, conformcount_raw_strict, col="slategray4", lwd=3.5, lty=1)
lines(xaxis, conformcount_raw_treat_strict, col="steelblue1", lwd=3.5, lty=5)
lines(xaxis, conformcount_raw_control_strict, col="black", lwd=3.5, lty=3)
axis(1, at=c(0,1,2,3,4,5,6,7,8,9,10), labels=xaxis)
legend(0,1, legend=c("All Participants", "Treatment", "Control"), lty=c(1, 5,3), lwd=rep(1,3), col=c("slategray4", "steelblue1", "black"), title="Cumulative Distribution", bty="n")
legend(6.5,.3, legend=c("All Participants", "Treatment", "Control"), lty=c(1,5,3), lwd=c(3.5, 2.5, 3.5), col=c("slategray4", "steelblue1", "black"), title="Raw Distribution", bty="n")

#Double check whether they're statistically distinguishable 
prop.test(x=c(conformcount_treat_strict2[10]*46, conformcount_control_strict2[10]*46), n=c(rep(length(pol2$conformcount_strict), 2)), p=NULL)

prop.test(x=c(conformcount_treat2[10]*63, conformcount_control2[10]*63), n=c(rep(length(pol$conformcount), 2)), p=NULL)

prop.test(x=c(conformcount_total_strict2[10]*46, conformcount_total2[10]*46), n=c(rep(length(pol2$conformcount_strict), 2)), p=NULL)

#Emotions during lab session (emotionplot)

happy<-table(pol$happy)
excited<-table(pol$excited)
angry<-table(pol$angry)
frustrated<-table(pol$frustrated)
anxious<-table(pol$anxiety)
confused<-table(pol$confused)
scared<-table(pol$scared)
surprised<-table(pol$surprised)

emotion<-c(happy, excited, angry, frustrated, anxious, confused, scared, surprised)
emotionprop<-c(emotion[1]/length(pol$happy), emotion[2]/length(pol$excited), emotion[3]/length(pol$angry), emotion[4]/length(pol$frustrated), emotion[5]/length(pol$anxiety), emotion[6]/length(pol$confused), emotion[7]/length(pol$scared), emotion[8]/length(pol$surprised))
emotionprop<-sort(emotionprop, decreasing=TRUE)
emotionprop

#surprise, anxiety, frustrated, confused, happy, excited, angry, scared

barplot(emotionprop, main="Emotional Experience in Lab", xlab="Emotion", ylab="Proportion of Participants", names.arg=c("Surprised", "Anxious", "Frustrated", "Confused", "Happy", "Excited", "Angry", "Scared"), col=c("lightskyblue", "lightskyblue", "lightskyblue", "lightskyblue", "lightskyblue", "lightskyblue", "lightskyblue", "lightskyblue"), ylim=c(0,0.50))

##Pressure Variables (sourcesplot)

#Plot where pressure is reported to come from
table(pol$pressure)
	#29 people reported having felt pressure
	29/46 #63% of those 46 who answered the question reported feeling the pressure
	t.test(pol$pressure ~ pol$condition) #no difference between treatment groups

friends<-table(pol$pressure_friends)
family<-table(pol$pressure_family)
boss<-table(pol$pressure_boss)
teacher<-table(pol$pressure_teacher)
classmates<-table(pol$pressure_classmates)

pressure<-c(friends/29, classmates/29, family/29, teacher/29, boss/29)
pressure
	#72% friends, 65.5% classmates, 62% family, 31% teacher, 13.8% boss
barplot(pressure, main="Sources of Political Pressure", xlab="Source", ylab="Proportion of Participants", names.arg=c("Friends", "Classmates", "Family", "Teacher", "Boss"), col=c("lightskyblue", "lightskyblue", "lightskyblue", "lightskyblue", "lightskyblue"), ylim=c(0,1.0))


#--------Balance Tables (in appendix)---------#
#Party ID
pidtest<-prop.test(
	x=c(length(pol$dem[pol$dem==1 & pol$condition==0]), length(pol$dem[pol$dem==1 & pol$condition==1])),
	n=c(length(pol$condition[pol$condition==0]), length(pol$condition[pol$condition==1])), 
	p=NULL)
pidrow<-c(pidtest$estimate, pidtest$p.value)

#Female
femaletest<-prop.test(
	x=c(length(pol$female[pol$female==1 & pol$condition==0]), length(pol$female[pol$female==1 & pol$condition==1])),
	n=c(length(pol$condition[pol$condition==0]), length(pol$condition[pol$condition==1])), 
	p=NULL)
femalerow<-c(femaletest$estimate, femaletest$p.value)
femalerow

#Race
whitetest<-prop.test(
	x=c(length(pol$white[pol$white==1 & pol$condition==0]), length(pol$white[pol$white==1 & pol$condition==1])),
	n=c(length(pol$condition[pol$condition==0]), length(pol$condition[pol$condition==1])), 
	p=NULL)
whiterow<-c(whitetest$estimate, whitetest$p.value)
whiterow


#2012 Engagement
engage2012test<-t.test(pol$engage2012 ~ pol$condition)
engage2012row<-c(engage2012test$estimate, engage2012test$p.value)
engage2012row

#Past 4 years engagement
engage4test<-t.test(pol$engage4 ~ pol$condition)
engage4row<-c(engage4test$estimate, engage4test$p.value)
engage4row

#Knowledge
knowtest<-t.test(pol$polknow ~ pol$condition)
knowrow<-c(knowtest$estimate, knowtest$p.value)
knowrow
table(pol$engage2012)
#Interest
interesttest<-t.test(pol$polinterest ~ pol$condition)
interestrow<-c(interesttest$estimate, interesttest$p.value)
interestrow

#Pay Attention
payattntest<-t.test(pol$payattn ~ pol$condition)
payattnrow<-c(payattntest$estimate, payattntest$p.value)
payattnrow

#Ideology (exclude 8 because that's haven't thought about it much)
ideologytest<-t.test(pol$ideology[pol$ideology!=8] ~ pol$condition)
ideologyrow<-c(ideologytest$estimate, ideologytest$p.value)
ideologyrow

#N
n<-c(length(pol$condition[pol$condition==0]), length(pol$condition[pol$condition==1]), 0)

study2balance<-matrix(NA, nrow=10, ncol=3)
study2balance[1,]<-pidrow
study2balance[2,]<-femalerow
study2balance[3,]<-whiterow
study2balance[4,]<-engage2012row
study2balance[5,]<-engage4row
study2balance[6,]<-knowrow
study2balance[7,]<-interestrow
study2balance[8,]<-payattnrow
study2balance[9,]<-ideologyrow
study2balance[10,]<-n
rownames(study2balance)<-c("Democrat", "Female", "White", "2012 Engagement", "Past 4 Years Engagement", "Knowledge", "Interest", "Pay Attention", "Ideology", "N")
colnames(study2balance)<-c("Control", "Treatment", "P-Value")

study2balance<-round(study2balance, digits=3)
library(xtable)
xtable(study2balance)

#---------Balance Table for Just Those Who Completed Posttest-------#

#Party ID
pidtest<-prop.test(
	x=c(length(pol2$dem[pol2$dem==1 & pol2$condition==0]), length(pol2$dem[pol2$dem==1 & pol2$condition==1])),
	n=c(length(pol2$condition[pol2$condition==0]), length(pol2$condition[pol2$condition==1])), 
	p=NULL)
pidrow<-c(pidtest$estimate, pidtest$p.value)

#Female
femaletest<-prop.test(
	x=c(length(pol2$female[pol2$female==1 & pol2$condition==0]), length(pol2$female[pol2$female==1 & pol2$condition==1])),
	n=c(length(pol2$condition[pol2$condition==0]), length(pol2$condition[pol2$condition==1])), 
	p=NULL)
femalerow<-c(femaletest$estimate, femaletest$p.value)
femalerow

#Race
whitetest<-prop.test(
	x=c(length(pol2$white[pol2$white==1 & pol2$condition==0]), length(pol2$white[pol2$white==1 & pol2$condition==1])),
	n=c(length(pol2$condition[pol2$condition==0]), length(pol2$condition[pol2$condition==1])), 
	p=NULL)
whiterow<-c(whitetest$estimate, whitetest$p.value)
whiterow


#2012 Engagement
engage2012test<-t.test(pol2$engage2012 ~ pol2$condition)
engage2012row<-c(engage2012test$estimate, engage2012test$p.value)
engage2012row

#Past 4 years engagement
engage4test<-t.test(pol2$engage4 ~ pol2$condition)
engage4row<-c(engage4test$estimate, engage4test$p.value)
engage4row

#Knowledge
knowtest<-t.test(pol2$polknow ~ pol2$condition)
knowrow<-c(knowtest$estimate, knowtest$p.value)
knowrow
table(pol2$engage2012)
#Interest
interesttest<-t.test(pol2$polinterest ~ pol2$condition)
interestrow<-c(interesttest$estimate, interesttest$p.value)
interestrow

#Pay Attention
payattntest<-t.test(pol2$payattn ~ pol2$condition)
payattnrow<-c(payattntest$estimate, payattntest$p.value)
payattnrow

#Ideology (exclude 8 because that's haven't thought about it much)
ideologytest<-t.test(pol2$ideology[pol2$ideology!=8] ~ pol2$condition)
ideologyrow<-c(ideologytest$estimate, ideologytest$p.value)
ideologyrow

#N
n<-c(length(pol2$condition[pol2$condition==0]), length(pol2$condition[pol2$condition==1]), 0)

study2balance<-matrix(NA, nrow=10, ncol=3)
study2balance[1,]<-pidrow
study2balance[2,]<-femalerow
study2balance[3,]<-whiterow
study2balance[4,]<-engage2012row
study2balance[5,]<-engage4row
study2balance[6,]<-knowrow
study2balance[7,]<-interestrow
study2balance[8,]<-payattnrow
study2balance[9,]<-ideologyrow
study2balance[10,]<-n
rownames(study2balance)<-c("Democrat", "Female", "White", "2012 Engagement", "Past 4 Years Engagement", "Knowledge", "Interest", "Pay Attention", "Ideology", "N")
colnames(study2balance)<-c("Control", "Treatment", "P-Value")

study2balance<-round(study2balance, digits=3)
library(xtable)
xtable(study2balance)

#--------------------------#


#Table comparing those who completed the full test to those who dropped out (in appendix)

#subset of the data of those who did not complete the posttest
nopost<-subset(pol, is.na(pol$pressure))
dim(nopost)

#Make a subset of nopost for a summary table (did not complete posttest)
nopost.sum<-nopost[, c(7:15, 17)]
library(stargazer)
stargazer(nopost.sum)

#Make a subset of pol2 for a summary table (completed posttest)
pol2.sum<-pol2[, c(7:15, 17)]
stargazer(pol2.sum)

#Make a subset of pol for a summary table (whole sample)
pol.sum<-pol[, c(7:15, 17)]
stargazer(pol.sum)



#---------Testing Hypotheses in Regression Framework---------#

#Run linear models

#Base model (potential)
pot.base<-lm(conformcount ~ condition, data=pol)
summary(pot.base)

#study artifacts (potential)
pot.study<-lm(conformcount ~ condition + remove + order, data=pol)
summary(pot.study)

#demographic controls (potential)
pot.demog<-lm(conformcount ~ condition + remove + order + female + white, data=pol)
summary(pot.demog)

#political controls (potential)
pot.pol<-lm(conformcount ~ condition + remove + order + female + white + engage2012 + polknow + polinterest + attach + ideology, data=pol)
summary(pot.pol)

#library(stargazer)
stargazer(pot.base, pot.study, pot.demog, pot.pol)

#Base model (pure)
pure.base<-lm(conformcount_strict ~ condition, data=pol2)
summary(pure.base)

#study artifacts (pure)
pure.study<-lm(conformcount_strict ~ condition + remove + order, data=pol2)
summary(pure.study)

#demographic controls (pure)
pure.demog<-lm(conformcount_strict ~ condition + remove + order + female + white, data=pol2)
summary(pure.demog)

#political controls (pure)
pure.pol<-lm(conformcount_strict ~ condition + remove + order + female + white + engage2012 + polknow + polinterest + attach + ideology, data=pol2)
summary(pure.pol)

stargazer(pure.base, pure.study, pure.demog, pure.pol)

#Poisson Model (in appendix)

#Base model (potential)
pot.base.pois<-glm(conformcount ~ condition, data=pol, family=poisson(link="log"))
summary(pot.base.pois)

#study artifacts (potential)
pot.study.pois<-glm(conformcount ~ condition + remove + order, data=pol, family=poisson(link="log"))
summary(pot.study.pois)

#demographic controls (potential)
pot.demog.pois<-glm(conformcount ~ condition + remove + order + female + white, data=pol, family=poisson(link="log"))
summary(pot.demog.pois)

#political controls (potential)
pot.pol.pois<-glm(conformcount ~ condition + remove + order + female + white + engage2012 + polknow + polinterest + attach + ideology, data=pol, family=poisson(link="log"))
summary(pot.pol.pois)
#library(stargazer)
stargazer(pot.base.pois, pot.study.pois, pot.demog.pois, pot.pol.pois)

#Base model (pure)
pure.base.pois<-glm(conformcount_strict ~ condition, data=pol2, family=poisson(link="log"))
summary(pure.base.pois)

#study artifacts (pure)
pure.study.pois<-glm(conformcount_strict ~ condition + remove + order, data=pol2, family=poisson(link="log"))
summary(pure.study.pois)

#demographic controls (pure)
pure.demog.pois<-glm(conformcount_strict ~ condition + remove + order + female + white, data=pol2, family=poisson(link="log"))
summary(pure.demog.pois)

#political controls (pure)
pure.pol.pois<-glm(conformcount_strict ~ condition + remove + order + female + white + engage2012 + polknow + polinterest + attach + ideology, data=pol2, family=poisson(link="log"))
summary(pure.pol.pois)

stargazer(pure.base.pois, pure.study.pois, pure.demog.pois, pure.pol.pois)









